Compiled under R version 4.3.0 (2023-04-21)
WARNING: edit the working directory to your preferred folder.
This document details all analyses performed in R for the study:
For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
library(tidyverse)
library(nlme)
library(ape)
library(phytools)
library(bioacoustics)
library(AICcmodavg)
library(evobiR)
library(phylopath)
library(caper)
library(geiger)
library(windex)
library(geomorph)
library(RColorBrewer)
fitEvolPar to estimate the best fit of
the main parameter (alpha, lambda, or g) in the corresponding
evolutionary model (OU, lambda, or EB, respectively)(Download from GitHub: https://github.com/LucasLegendre/fitEvolPar)
sys.source("fitEvolPar.R", envir = knitr::knit_global())
# Load both datasets
acoudata<-read.csv("notes-filtered_20211202.csv", header=TRUE, sep=";")
datasyr<-read.table("Dataset_syrinx_strisores.txt", header=TRUE)
rownames(datasyr)<-datasyr$Taxon
# Extract species in acoustic dataset for which we have morphological data
acoudataset<-acoudata %>%
group_by(species) %>%
dplyr::summarize(duration=mean(duration), freq_max_amp=mean(freq_max_amp), freq_min=mean(freq_min),
freq_max=mean(freq_max), bandwidth=mean(bandwidth))
acoudataset<-as.data.frame(acoudataset); rownames(acoudataset)<-acoudataset$species
vector<-c(datasyr$Taxon[1], "Chaetura_pelagica",
datasyr$Taxon[c(4:12)],"Oreotrochilus_chimborazo", datasyr$Taxon[c(15:21)])
intersect(acoudataset$species, vector) # Species present in both datasets
## [1] "Apus_affinis" "Archilochus_colubris"
## [3] "Calliphlox_amethystina" "Chaetura_pelagica"
## [5] "Chordeiles_minor" "Coeligena_violifer"
## [7] "Collocalia_esculenta" "Eutoxeres_condamini"
## [9] "Florisuga_mellivora" "Heliangelus_amethysticollis"
## [11] "Mellisuga_minima" "Oreotrochilus_chimborazo"
## [13] "Patagona_gigas" "Phaethornis_guy"
## [15] "Phaethornis_hispidus" "Phaethornis_superciliosus"
## [17] "Topaza_pella"
acoumorph<-acoudataset %>% filter(species %in% vector)
# Reformat to match syrinx data
acoumorph<-acoumorph[match(vector,acoumorph$species),] # reorder to match syrinx data
rownames(acoumorph)[c(17,18)]<-datasyr$Taxon[c(19,20)]
acoumorph$species[c(17,18)]<-rownames(acoumorph)[c(17,18)]
acoumorph<-acoumorph[c(1,rep(2,2),3:11,rep(12,2),13:19),]
# Bind two datasets
data<-cbind(datasyr,acoumorph[,c(2:6)])
data<-data[,c(1:11,22:27)]; colnames(data)[12]<-"IM_CSA"
rownames(data)<-data$Taxon
Here, we extract five acoustic traits to add them to the
morphological dataset – duration, frequency at maximum amplitude,
minimum frequency, maximum frequency, and bandwidth – for more info on
each parameter, see vignette for package bioacoustics:
vignette("introduction", package = "bioacoustics")
## starting httpd help server ... done
treeS<-read.nexus("Tree_strisores.trees.nex")
treesyr<-drop.tip(treeS, setdiff(treeS$tip.label, rownames(datasyr)))
plotTree(treesyr, fsize=1.2)
data<-ReorderData(treesyr,data)
# Remove missing data
acouphyl<-data[!is.na(data$duration),c(13:17)]
acoutree<-drop.tip(treesyr, setdiff(treesyr$tip.label,rownames(acouphyl)))
# Phylogenetic signal
var=list(); phy=list()
for (i in 1:5) {
var<-acouphyl[,i]; names(var)<-rownames(acouphyl)
phy[[i]]<-phylosig(acoutree, var, method="lambda", test=TRUE)
}
phy
## [[1]]
##
## Phylogenetic signal lambda : 1.00156
## logL(lambda) : -87.9944
## LR(lambda=0) : 17.3421
## P-value (based on LR test) : 3.12193e-05
##
##
## [[2]]
##
## Phylogenetic signal lambda : 6.62145e-05
## logL(lambda) : -166.581
## LR(lambda=0) : -0.000438143
## P-value (based on LR test) : 1
##
##
## [[3]]
##
## Phylogenetic signal lambda : 6.62145e-05
## logL(lambda) : -161.255
## LR(lambda=0) : -0.000812817
## P-value (based on LR test) : 1
##
##
## [[4]]
##
## Phylogenetic signal lambda : 1.00156
## logL(lambda) : -166.028
## LR(lambda=0) : 13.6852
## P-value (based on LR test) : 0.00021615
##
##
## [[5]]
##
## Phylogenetic signal lambda : 1.00156
## logL(lambda) : -158.831
## LR(lambda=0) : 13.1823
## P-value (based on LR test) : 0.000282605
The signal is strong for three variables (duration, max. frequency, and bandwidth), but not significant for the other two (frequency at max. amplitude and min. frequency). We will use a lambda model (Brownian Motion with Pagel’s lambda optimized using restricted maximum likelihood).
phytools)pPCA<-phyl.pca(acoutree, acouphyl, method="lambda", opt="REML")
## [1] 0.3825872 -234.7789361
## [1] 0.6190391 -245.5371293
## [1] 0.2364519 -229.5685693
## [1] 0.1461353 -227.2380061
## [1] 0.09031659 -225.74346049
## [1] 0.05581872 -224.85254229
## [1] 0.03449787 -224.36620861
## [1] 0.02132085 -223.95581601
## [1] 0.01317701 -224.00769067
## [1] 0.01905903 -223.86355678
## [1] 0.01764617 -224.31072533
## [1] 0.01992297 -224.11316405
## [1] 0.01851937 -224.03902107
## [1] 0.01938903 -224.03286074
## [1] 0.0188529 -223.9688029
## [1] 0.01918508 -224.09870716
## [1] 0.0189803 -223.9723786
## [1] 0.01910718 -224.14145637
## [1] 0.01905903 -223.86355678
summary(pPCA)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 401.2734513 141.2691511 45.96595318 4.6742621244
## Proportion of Variance 0.8793543 0.1089878 0.01153866 0.0001193189
## Cumulative Proportion 0.8793543 0.9883420 0.99988068 1.0000000000
## PC5
## Standard deviation 1.749780e-05
## Proportion of Variance 1.672053e-15
## Cumulative Proportion 1.000000e+00
plot(pPCA); biplot(pPCA)
pPCA$lambda # very low value, as expected
## [1] 0.01905903
PC1 explains 88% of the variance, PC2 only 11% – we can use PC1 as a character representing acoustic parameters for the species in our dataset.
PC1<-pPCA$S[,1]
PC1<-c(PC1[1:12],NA,PC1[13:17],NA,PC1[18:19])
names(PC1)[c(13,19)]<-rownames(acoumorph)[c(20,19)]
data<-cbind(data, PC1)
var=list(); phy=list(); phydat=data.frame()
for (i in 4:12) {
var<-data[,i]; names(var)<-rownames(data)
var2<-var[!is.na(var)]
treesyrvar<-drop.tip(treesyr,setdiff(treesyr$tip.label,names(var2)))
phy[[i]]<-phylosig(treesyrvar, var2, method="lambda", test=TRUE)
phydat[i,1]<-colnames(data)[i]; phydat[i,2]<-phy[[i]]$lambda; phydat[i,3]<-phy[[i]]$P
}
colnames(phydat)<-c("Variable", "Lambda", "P")
phydat[c(4:12),]
### Checking evolutionary model with fitContinuous
models=c("BM", "OU", "EB", "rate_trend","lambda", "white") # (For more info, check '?fitContinuous')
var=list(); fit=list(); mod=list()
for (i in 4:12) {
var<-log(data[,i]); names(var)<-rownames(data)
var2<-var[!is.na(var)]
treesyrvar<-drop.tip(treesyr,setdiff(treesyr$tip.label,names(var2)))
for (m in 1:length(models)) {
fit[[m]]=fitContinuous(treesyrvar, var2, model=models[m], ncores=2)
}
mod[[i]]<-modSel.geiger(fit[[1]],fit[[2]],fit[[3]],fit[[4]],fit[[5]],fit[[6]])
}
mod
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[4]] 3 49.17265 -90.93353 0.0000 1 1.000000e+00
## fit[[1]] 2 -15.00120 34.66906 125.6026 0 1.880423e+27
## fit[[3]] 3 -14.94029 37.29234 128.2259 0 6.980620e+27
## fit[[2]] 3 -15.00120 37.41416 128.3477 0 7.419012e+27
## fit[[5]] 3 -15.00120 37.41416 128.3477 0 7.419012e+27
## fit[[6]] 2 -27.37592 59.41850 150.3520 0 4.451727e+32
##
## [[5]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[6]] 2 -4.082598 12.87108 0.000000 0.5580 1.000000e+00
## fit[[5]] 3 -3.539262 14.57852 1.707445 0.2376 2.348373e+00
## fit[[2]] 3 -3.689804 14.87961 2.008529 0.2044 2.729899e+00
## fit[[4]] 3 -13.985357 35.47071 22.599637 0.0000 8.080696e+04
## fit[[1]] 2 -15.925669 36.55722 23.686142 0.0000 1.391170e+05
## fit[[3]] 3 -15.925762 39.35152 26.480445 0.0000 5.625428e+05
##
## [[6]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[2]] 3 -16.94350 41.29877 0.000000 0.7127 1.000000
## fit[[5]] 3 -18.63204 44.67584 3.377070 0.1317 5.411546
## fit[[6]] 2 -20.71333 46.09333 4.794562 0.0648 10.993246
## fit[[4]] 3 -19.78713 46.98603 5.687262 0.0415 17.178030
## fit[[1]] 2 -21.21266 47.09199 5.793215 0.0393 18.112595
## fit[[3]] 3 -21.21272 49.83721 8.538440 0.0100 71.465871
##
## [[7]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[5]] 3 -12.65971 32.73118 0.000000 0.7849 1.000000
## fit[[2]] 3 -14.63783 36.68742 3.956238 0.1086 7.229132
## fit[[1]] 2 -16.84765 38.36196 5.630776 0.0470 16.699653
## fit[[4]] 3 -15.72212 38.85600 6.124818 0.0367 21.378999
## fit[[3]] 3 -16.84770 41.10717 8.375988 0.0119 65.890493
## fit[[6]] 2 -18.31153 41.28973 8.558546 0.0109 72.187955
##
## [[8]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[1]] 2 -9.212933 23.28301 0.000000 0.3446 1.000000
## fit[[2]] 3 -8.366069 24.57829 1.295283 0.1803 1.911028
## fit[[4]] 3 -8.572344 24.99084 1.707832 0.1467 2.348827
## fit[[5]] 3 -8.650226 25.14661 1.863597 0.1357 2.539071
## fit[[6]] 2 -10.306662 25.47047 2.187458 0.1154 2.985386
## fit[[3]] 3 -9.212967 26.27209 2.989079 0.0773 4.457284
##
## [[9]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[1]] 2 -8.041426 20.94000 0.000000 0.5265 1.000000
## fit[[4]] 3 -8.040869 23.92789 2.987897 0.1182 4.454649
## fit[[3]] 3 -8.040938 23.92803 2.988034 0.1182 4.454956
## fit[[2]] 3 -8.041426 23.92901 2.989011 0.1181 4.457132
## fit[[5]] 3 -8.041426 23.92901 2.989011 0.1181 4.457132
## fit[[6]] 2 -14.441956 33.74105 12.801059 0.0009 602.163690
##
## [[10]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[1]] 2 1.750625 1.204632 0.0000000 0.3603 1.000000
## fit[[2]] 3 2.826509 1.846981 0.6423492 0.2613 1.378746
## fit[[4]] 3 2.509289 2.481421 1.2767893 0.1903 1.893439
## fit[[5]] 3 1.750625 3.998750 2.7941176 0.0891 4.043290
## fit[[3]] 3 1.750591 3.998817 2.7941854 0.0891 4.043428
## fit[[6]] 2 -1.845972 8.397826 7.1931945 0.0099 36.473912
##
## [[11]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[2]] 3 1.2034974 5.093005 0.000000 0.4416 1.000000
## fit[[1]] 2 -0.7733693 6.252621 1.159616 0.2473 1.785695
## fit[[4]] 3 0.2556109 6.988778 1.895773 0.1711 2.580250
## fit[[5]] 3 -0.7016153 8.903231 3.810225 0.0657 6.720164
## fit[[3]] 3 -0.7734162 9.046832 3.953827 0.0612 7.220423
## fit[[6]] 2 -3.7121227 12.130128 7.037122 0.0131 33.735856
##
## [[12]]
## K logLik AICc deltaAICc Weight Evidence ratio
## fit[[6]] 2 -14.60581 34.21162 0.000000 0.7097 1.000000e+00
## fit[[5]] 3 -14.44931 37.08043 2.868815 0.1691 4.197157e+00
## fit[[2]] 3 -14.78187 37.74556 3.533948 0.1212 5.853116e+00
## fit[[4]] 3 -24.22837 56.63856 22.426944 0.0000 7.412232e+04
## fit[[1]] 2 -27.30210 59.60419 25.392577 0.0000 3.265336e+05
## fit[[3]] 3 -27.30215 62.78611 28.574495 0.0000 1.602776e+06
High and significant values of lambda for all traits except distance TL-labia and IM CSA.
Using phylopath.
# Trachea length as % trachea length/total vocal tract length
TracheaRatio<-(data[,8]/(data[,8]+data[,9]))*100; names(TracheaRatio)<-rownames(data)
PPAdata<-cbind(data[,c(1,4:7)], TracheaRatio, data[,c(10,11,18)])
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
# For more information, see '?corClasses' in ape
# With body mass as predictor
bodymass<-PPAdata$Body_mass; names(bodymass)<-paste(rownames(PPAdata))
for (i in (1:7)) {
var2<-log10(PPAdata[,i+2])
names(var2)<-rownames(PPAdata)
var<-subset(var2, !is.na(var2))
treetest<-drop.tip(treesyr, setdiff(treesyr$tip.label, names(var)))
bodymass2<-log10(subset(bodymass, !is.na(var2)))
datatest<-as.data.frame(cbind(bodymass2,var))
spp<-names(var)
model=list()
model[[1]]=gls(var~bodymass2, datatest,
correlation=corBrownian(phy=treetest, form=~spp), method="ML")
model[[2]]=gls(var~bodymass2, datatest,
correlation=corMartins(fitEvolPar(datatest, treetest,"OU"),
phy=treetest, fixed=TRUE,
form=~spp), method="ML")
model[[3]]=gls(var~bodymass2, datatest,
correlation=corPagel(fitEvolPar(datatest,treetest,"lambda"),
phy=treetest,
fixed=TRUE, form=~spp), method="ML")
model[[4]]=gls(var~bodymass2, datatest,
correlation=corBlomberg(fitEvolPar(datatest,treetest,"EB"),
phy=treetest, fixed=TRUE,
form=~spp), method="ML")
model[[5]]=gls(var~bodymass2, datatest, method="ML")
print(colnames(PPAdata)[i+2])
print(aictab(cand.set=model, modnames=Modnames,sort=TRUE))
}
## [1] 0
## [1] 0.6
## [1] 0.1
## [1] "Distance_TL.labia"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 5.69 Inf NaN NaN 0.90
## Lambda 3 -20.62 Inf NaN NaN 14.06
## EB 3 -9.86 Inf NaN NaN 8.68
## OLS 3 -17.70 Inf NaN NaN 12.60
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0.9
## [1] 0.1
## [1] "TL_CSA"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 5.56 Inf NaN NaN 0.93
## Lambda 3 -7.52 Inf NaN NaN 7.46
## EB 3 -2.15 Inf NaN NaN 4.78
## OLS 3 1.74 Inf NaN NaN 2.84
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0.9
## [1] 0.2
## [1] "Length_labia"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 3.57 Inf NaN NaN 1.92
## Lambda 3 -3.54 Inf NaN NaN 5.47
## EB 3 1.79 Inf NaN NaN 2.81
## OLS 3 8.70 Inf NaN NaN -0.64
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0.9
## [1] 1
## [1] "TracheaRatio"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 -28.09 Inf NaN NaN 17.97
## Lambda 3 -29.11 Inf NaN NaN 18.48
## EB 3 -28.09 Inf NaN NaN 17.97
## OLS 3 -20.54 Inf NaN NaN 14.19
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "Tracheal_diameter_EXT"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 -41.79 Inf NaN NaN 24.65
## Lambda 3 -54.22 Inf NaN NaN 30.86
## EB 3 -53.07 Inf NaN NaN 30.29
## OLS 3 -54.22 Inf NaN NaN 30.86
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "Tracheal_diameter_INT"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 -32.74 Inf NaN NaN 20.12
## Lambda 3 -48.83 Inf NaN NaN 28.17
## EB 3 -46.65 Inf NaN NaN 27.07
## OLS 3 -48.83 Inf NaN NaN 28.17
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "PC1"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 22.32 Inf NaN NaN -6.16
## Lambda 3 12.24 Inf NaN NaN -1.12
## EB 3 12.05 Inf NaN NaN -1.03
## OLS 3 12.24 Inf NaN NaN -1.12
## OU 3 -Inf NaN NaN NaN Inf
# With PC1 as response
for (i in (1:7)) {
var2<-log10(PPAdata[,i+1])
names(var2)<-rownames(PPAdata)
var<-subset(var2, !is.na(var2)&!is.na(PC1))
PC12<-subset(PC1, !is.na(PC1)&!is.na(var2))
treetest<-drop.tip(treesyr, setdiff(treesyr$tip.label, names(var)))
datatest<-as.data.frame(cbind(var,PC12))
spp<-names(var)
model=list()
model[[1]]=gls(PC12~var, datatest,
correlation=corBrownian(phy=treetest, form=~spp), method="ML")
model[[2]]=gls(PC12~var, datatest,
correlation=corMartins(fitEvolPar(datatest, treetest,"OU"),
phy=treetest, fixed=TRUE,
form=~spp), method="ML")
model[[3]]=gls(PC12~var, datatest,
correlation=corPagel(fitEvolPar(datatest,treetest,"lambda"),
phy=treetest,
fixed=TRUE, form=~spp), method="ML")
model[[4]]=gls(PC12~var, datatest,
correlation=corBlomberg(fitEvolPar(datatest,treetest,"EB"),
phy=treetest, fixed=TRUE,
form=~spp), method="ML")
model[[5]]=gls(PC12~var, datatest, method="ML")
print(colnames(PPAdata)[i+1])
print(aictab(cand.set=model, modnames=Modnames,sort=TRUE))
}
## [1] 0
## [1] 1
## [1] 0.1
## [1] "Body_mass"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 358.02 Inf NaN NaN -175.21
## Lambda 3 358.02 Inf NaN NaN -175.21
## EB 3 352.54 Inf NaN NaN -172.47
## OLS 3 364.98 Inf NaN NaN -178.69
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 1
## [1] 0.1
## [1] "Distance_TL.labia"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 340.23 Inf NaN NaN -166.26
## Lambda 3 340.23 Inf NaN NaN -166.26
## EB 3 335.82 Inf NaN NaN -164.05
## OLS 3 347.16 Inf NaN NaN -169.72
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 1
## [1] 0.1
## [1] "TL_CSA"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 358.76 Inf NaN NaN -175.58
## Lambda 3 358.76 Inf NaN NaN -175.58
## EB 3 354.47 Inf NaN NaN -173.43
## OLS 3 366.41 Inf NaN NaN -179.41
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 1
## [1] 0.1
## [1] "Length_labia"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 358.88 Inf NaN NaN -175.64
## Lambda 3 358.88 Inf NaN NaN -175.64
## EB 3 354.71 Inf NaN NaN -173.56
## OLS 3 366.37 Inf NaN NaN -179.39
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "TracheaRatio"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 288.38 Inf NaN NaN -140.10
## Lambda 3 284.45 Inf NaN NaN -138.13
## EB 3 284.50 Inf NaN NaN -138.16
## OLS 3 284.45 Inf NaN NaN -138.13
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "Tracheal_diameter_EXT"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 345.05 Inf NaN NaN -168.67
## Lambda 3 341.82 Inf NaN NaN -167.05
## EB 3 338.86 Inf NaN NaN -165.57
## OLS 3 341.82 Inf NaN NaN -167.05
## OU 3 -Inf NaN NaN NaN Inf
##
## [1] 0
## [1] 0
## [1] 0.1
## [1] "Tracheal_diameter_INT"
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## BM 3 345.56 Inf NaN NaN -168.93
## Lambda 3 342.37 Inf NaN NaN -167.33
## EB 3 339.50 Inf NaN NaN -165.89
## OLS 3 342.37 Inf NaN NaN -167.33
## OU 3 -Inf NaN NaN NaN Inf
In most cases, the BM model is selected as the best fit – we will use it for PPA.
Due to the high number of variables (Hardenberg and Gonzales-Voyer, 2013; Gonzales-Voyer and Hardenberg, 2014), we will decompose the PPA into two separate analyses: one on the effect of body mass on syringeal variables, and one on the effect of syringeal variables on sound production (represented by PC1 from our pPCA).
We can remove one of the two tracheal diameters, since the ratio of the two remains fairly constant through our sample:
RatioTD<-(log10(PPAdata$Tracheal_diameter_INT)/log10(PPAdata$Tracheal_diameter_EXT))*100
RatioTD<-RatioTD[1:20]
min(RatioTD); max(RatioTD) # values between 95 and 99% when log-converted
## [1] 95.01686
## [1] 98.70606
We will be using the external diameter in subsequent analyses.
PPAdata[,c(2:8)]<-log10(PPAdata[,c(2:8)])
colnames(PPAdata)[2:8]<-c("BM","DTL","CSA","LL","Tr","TDe","TDi")
M<-define_model_set(null=c(),
one=c(LL~Tr),
two=c(LL~Tr, Tr~BM),
three=c(LL~TDe),
four=c(LL~TDe+Tr),
five=c(LL~TDe+Tr, Tr~BM),
six=c(LL~DTL),
seven=c(LL~DTL, DTL~CSA),
eight=c(LL~DTL, DTL~CSA, CSA~Tr),
nine=c(LL~DTL, DTL~CSA, CSA~Tr, Tr~BM),
ten=c(LL~DTL, DTL~CSA, CSA~TDe),
eleven=c(LL~DTL, DTL~CSA, CSA~TDe+Tr),
twelve=c(LL~DTL, DTL~CSA, CSA~TDe+Tr, Tr~BM),
thirteen=c(LL~DTL+Tr),
fourteen=c(LL~DTL+Tr, DTL~CSA),
fifteen=c(LL~DTL+Tr, DTL~CSA, CSA~Tr),
sixteen=c(LL~DTL+Tr, DTL~CSA, CSA~Tr, Tr~BM),
seventeen=c(LL~DTL+Tr, DTL~CSA, CSA~TDe),
eighteen=c(LL~DTL+Tr, DTL~CSA, CSA~TDe+Tr),
nineteen=c(LL~DTL+Tr, DTL~CSA, CSA~TDe+Tr, Tr~BM),
twenty=c(LL~DTL+Tr, Tr~BM),
twentytwo=c(LL~DTL+Tr, DTL~CSA, Tr~BM),
twentythree=c(LL~DTL+Tr, DTL~CSA, CSA~Tr, Tr~BM),
twentyfour=c(LL~DTL+Tr, DTL~CSA, CSA~Tr, Tr~BM),
twentyfive=c(LL~DTL+Tr, DTL~CSA, CSA~TDe, Tr~BM),
twentysix=c(LL~DTL+Tr, DTL~CSA, CSA~TDe+Tr, Tr~BM),
twentyseven=c(LL~DTL+Tr, DTL~CSA, CSA~TDe+Tr, Tr~BM),
twentyeight=c(LL~DTL+TDe),
twentynine=c(LL~DTL+TDe, DTL~CSA),
thirty=c(LL~DTL+TDe, DTL~CSA, CSA~Tr, Tr~BM),
thirtyone=c(LL~DTL+TDe, DTL~CSA, CSA~TDe),
thirtytwo=c(LL~DTL+TDe, DTL~CSA, CSA~TDe+Tr),
thirtythree=c(LL~DTL+TDe, DTL~CSA, CSA~TDe+Tr, Tr~BM),
thirtyfour=c(LL~DTL+Tr+TDe),
thirtyfive=c(LL~DTL+Tr+TDe, DTL~CSA),
thirtysix=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~Tr),
thirtyseven=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~Tr, Tr~BM),
thirtyeight=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~TDe),
thirtynine=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~TDe+Tr),
forty=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~TDe+Tr, Tr~BM),
fortyone=c(LL~DTL+Tr+TDe, Tr~BM),
fortytwo=c(LL~DTL+Tr+TDe, DTL~CSA, Tr~BM),
fortythree=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~Tr, Tr~BM),
fortyfour=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~TDe, Tr~BM),
fortyfive=c(LL~DTL+Tr+TDe, DTL~CSA, CSA~TDe+Tr, Tr~BM),
.common=c(TDi~TDe,CSA~BM,DTL~BM,LL~BM,TDe~BM,TDi~BM, LL~BM))
plot_model_set(M)
result<-phylo_path(M, data=PPAdata, tree=treesyr, model='BM', "logistic_IG10"); result
## 5 rows were dropped because they contained NA values.
## Pruned tree to drop species not included in dat.
## A phylogenetic path analysis, on the variables:
## Continuous: Tr BM LL DTL CSA TDe TDi
## Binary:
##
## Evaluated for these models: null one two three four five six seven eight nine ten eleven twelve thirteen fourteen fifteen sixteen seventeen eighteen nineteen twenty twentytwo twentythree twentyfour twentyfive twentysix twentyseven twentyeight twentynine thirty thirtyone thirtytwo thirtythree thirtyfour thirtyfive thirtysix thirtyseven thirtyeight thirtynine forty fortyone fortytwo fortythree fortyfour fortyfive
##
## Containing 492 phylogenetic regressions, of which 53 unique
s<-summary(result); s
## Warning in summary.phylopath(result): CICc was not calculated for causal models
## where the number of parameters is equal toor larger than the number of species.
plot(s)
## Warning: Removed 45 rows containing missing values (`position_stack()`).
## Warning: Removed 45 rows containing missing values (`geom_text()`).
Not enough data to compile CICc due to small sample size – we will focus on individual regressions.
varpairs<-combn(c(colnames(PPAdata[2:7]),"IM_CSA","PC1"),2)
varpairs<-varpairs[,-c(11,14,23)] # Hypotheses we are not testing
varpair<-list()
for (i in c(1:7,9,11,12,15,16,20:25)) {
varpair[[i]]<-paste(varpairs[2,i], varpairs[1,i], sep="~")
}
for (i in c(8,10,13,14,17:19)) {
varpair[[i]]<-paste(varpairs[1,i], varpairs[2,i], sep="~")
}
# Create new dataset to include intrinsic muscle CSA
IM_CSA<-log10(data$IM_CSA)
datanew<-cbind(PPAdata, IM_CSA)
# PGLS
modlist<-as.data.frame(matrix(NA, nrow=length(varpair), ncol=3))
colnames(modlist)<-c("Model","R2","p")
for (i in c(1:19,21:25)) {
datapair<-cbind.data.frame(datanew$Taxon, datanew[,paste(varpairs[1,i])],
datanew[,paste(varpairs[2,i])])
colnames(datapair)<-c("Taxon",varpairs[1,i],varpairs[2,i])
comppair<-comparative.data(treesyr, datapair, names.col="Taxon")
modpair<-pgls(as.formula(paste(varpair[[i]])), data=comppair, lambda="ML")
modlist[i,1]<-paste(varpair[[i]])
modlist[i,2]<-summary(modpair)$adj.r.squared
modlist[i,3]<-summary(modpair)$coef[2,4]
}
# Model 20 (PC1~LL) cannot compute because of problem with 'optim' (ML estimate of lambda is probably out of bounds)
# We can get the value of lambda with 'fitEvolPar', compute the corresponding model, and add it to the data frame of models
data20<-cbind.data.frame(datanew[,paste(varpairs[1,20])],
datanew[,paste(varpairs[2,20])])
colnames(data20)<-c(varpairs[1,20],varpairs[2,20])
rownames(data20)<-rownames(PPAdata)
data20<-na.omit(data20)
tree20<-drop.tip(treesyr, setdiff(treesyr$tip.label, rownames(data20)))
fitEvolPar(data20, tree20, "lambda")
## [1] 1
# lambda = 1: we do not need an ML estimate for lambda
for (i in 20) {
datapair<-cbind.data.frame(datanew$Taxon, datanew[,paste(varpairs[1,i])],
datanew[,paste(varpairs[2,i])])
colnames(datapair)<-c("Taxon",varpairs[1,i],varpairs[2,i])
comppair<-comparative.data(treesyr, datapair, names.col="Taxon")
modpair<-pgls(as.formula(paste(varpair[[i]])), data=comppair)
modlist[i,1]<-paste(varpair[[i]])
modlist[i,2]<-summary(modpair)$adj.r.squared
modlist[i,3]<-summary(modpair)$coef[2,4]
}
modlist # list of all models with pseudo R-squared and p-values
Let us check which models are significant
modlist[which(modlist$p<0.05),]
for (i in which(modlist$p<0.05)) {
datapair<-cbind.data.frame(datanew$Taxon, datanew[,paste(varpairs[1,i])],
datanew[,paste(varpairs[2,i])])
colnames(datapair)<-c("Taxon",varpairs[1,i],varpairs[2,i])
treepair<-drop.tip(treesyr, setdiff(treesyr$tip.label, datapair$Taxon))
comppair<-comparative.data(treepair, datapair, names.col="Taxon")
modpair<-pgls(as.formula(paste(varpair[[i]])), data=comppair, lambda="ML")
print(as.formula(paste(varpair[[i]])))
print(summary(modpair)$coef)
print(paste("R-squared =", summary(modpair)$adj.r.squared))
}
## CSA ~ BM
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8809359 0.2579353 15.046163 5.214496e-12
## BM 0.7983534 0.1519696 5.253376 4.528547e-05
## [1] "R-squared = 0.570796661549721"
## LL ~ BM
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8421281 0.2776243 6.635327 2.392372e-06
## BM 0.3618102 0.1620229 2.233081 3.776862e-02
## [1] "R-squared = 0.16620292803003"
## TDe ~ BM
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9380373 0.03027184 97.055114 0.000000e+00
## BM 0.2597343 0.03073778 8.450002 1.114473e-07
## [1] "R-squared = 0.787478059507324"
## IM_CSA ~ BM
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5296253 0.1762632 31.371420 1.223466e-13
## BM 0.7066604 0.2148268 3.289442 5.865668e-03
## [1] "R-squared = 0.412269251188365"
## DTL ~ CSA
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9594961 0.4167421 4.701940 0.0001776936
## CSA 0.2658761 0.0855205 3.108917 0.0060599328
## [1] "R-squared = 0.313220671406667"
## DTL ~ Tr
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.749126 0.3929901 12.084594 8.546387e-09
## Tr -0.840497 0.2238408 -3.754888 2.132676e-03
## [1] "R-squared = 0.466176663866154"
## IM_CSA ~ DTL
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1241414 1.4212382 -0.08734738 0.9318360597
## DTL 1.8648043 0.4250838 4.38691027 0.0008854953
## [1] "R-squared = 0.583933186384203"
## CSA ~ TDe
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.826354 1.2971701 -1.407953 1.761818e-01
## TDe 2.071955 0.4014746 5.160862 6.563061e-05
## [1] "R-squared = 0.574320293768139"
## IM_CSA ~ CSA
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8405865 0.9615453 2.954189 0.011181491
## CSA 0.6722607 0.1993087 3.372962 0.004995249
## [1] "R-squared = 0.425685221282917"
## LL ~ Tr
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.931236 0.7850172 6.281692 1.470521e-05
## Tr -1.412626 0.4479796 -3.153327 6.561557e-03
## [1] "R-squared = 0.358549580541089"
## LL ~ IM_CSA
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8336393 0.37938540 2.197342 0.0467250578
## IM_CSA 0.2828725 0.06017951 4.700479 0.0004148586
## [1] "R-squared = 0.601077116895975"
## IM_CSA ~ TDe
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5083212 1.6185512 -0.3140594 0.758868645
## TDe 2.0997796 0.5135103 4.0890702 0.001501441
## [1] "R-squared = 0.547361569552698"
Tracheal length ratio is negatively correlated with distance TL-labia, which is positively correlated with CSA of intrinsic muscles. In hummingbirds, shortening of the trachea may have shortened TL and liberated extra space for intrinsic muscles to develop at the level of the tympanum. Additionally, both shortening of the trachea and higher CSA of intrinsic muscles are correlated with longer labia, which may have an effect on hummingbird vocal production. However, only one of our investigated traits, shortening of the trachea, is correlated with PC1 from our pPCA on sound characters. It is possible that the effect of tracheal length ratio on CSA of intrinsic muscles has an effect on sound production undetected by our analyses due to small sample size; this hypothesis, however, requires further testing.
We can test this hypothesis by including individual acoustic traits and correlating them with morphological traits.
dataS<-cbind(datanew[,c(1:7,10)],log10(data[,c(13:17)]))
# List pairs of variables (x, y) to compile regressions of
pairsound<-expand.grid(colnames(dataS)[2:8],colnames(dataS)[9:13])
# Perform all PGLS
modlistsound<-as.data.frame(matrix(NA, nrow=35, ncol=3))
colnames(modlistsound)<-c("Model","R2","p")
for (i in c(1:10,12:35)) {
datapgls<-cbind.data.frame(dataS$Taxon, dataS[,paste(pairsound[i,1])],
dataS[,paste(pairsound[i,2])])
colnames(datapgls)<-c("Taxon",paste(pairsound[i,1]),paste(pairsound[i,2]))
compS<-comparative.data(treesyr, datapgls, names.col="Taxon")
pglS<-pgls(as.formula(paste(colnames(datapgls[3]),"~",colnames(datapgls[2]))),
data=compS, lambda="ML")
modlistsound[i,1]<-as.character(paste(colnames(datapgls[3]),"~",colnames(datapgls[2])))
modlistsound[i,2]<-summary(pglS)$adj.r.squared
modlistsound[i,3]<-summary(pglS)$coef[2,4]
}
# Model 20 (freq_max_amp~LL) cannot compute because of problem with 'optim' (ML estimate of lambda is probably out of bounds)
# We can get the value of lambda with 'fitEvolPar', compute the corresponding model, and add it to the data frame of models
i = 11
datapgls<-cbind.data.frame(dataS[,paste(pairsound[i,1])],
dataS[,paste(pairsound[i,2])])
colnames(datapgls)<-c(paste(pairsound[i,1]),paste(pairsound[i,2]))
rownames(datapgls)<-dataS$Taxon
datapgls<-na.omit(datapgls)
tree11<-drop.tip(treesyr, setdiff(treesyr$tip.label, rownames(datapgls)))
fitEvolPar(datapgls, tree11, "lambda")
## [1] 1
# lambda = 1: we do not need an ML estimate for lambda
datapgls<-cbind.data.frame(dataS$Taxon,dataS[,paste(pairsound[i,1])],
dataS[,paste(pairsound[i,2])])
colnames(datapgls)<-c("Taxon", paste(pairsound[i,1]),paste(pairsound[i,2]))
compS<-comparative.data(treesyr, datapgls, names.col="Taxon")
pglS<-pgls(as.formula(paste(colnames(datapgls[3]),"~",colnames(datapgls[2]))),
data=compS)
modlistsound[i,1]<-as.character(paste(colnames(datapgls[3]),"~",colnames(datapgls[2])))
modlistsound[i,2]<-summary(pglS)$adj.r.squared
modlistsound[i,3]<-summary(pglS)$coef[2,4]
# Check for significance
min(modlistsound$R2); max(modlistsound$R2)
## [1] -0.08708794
## [1] 0.1983234
which(modlistsound$p<0.05)
## integer(0)
None of the regressions are significant, probably due to small sample size.
# Data matrices
dataS2<-dataS[,c(2:7,9:13)] # No IM CSA due to a high amount of missing data (all non-hummingbirds)
dataS2<-na.omit(dataS2)
Morphol<-dataS2[,1:6]
Sound<-dataS2[,7:11]
treepls<-drop.tip(treesyr, setdiff(treesyr$tip.label, rownames(dataS2)))
# 2B-PLS
twobpls<-phylo.integration(log1p(Morphol), log1p(Sound), phy=treepls, iter=9999) # with phylogeny
summary(twobpls); plot(twobpls)
twobpls2<-two.b.pls(log1p(Morphol), log1p(Sound), iter=9999) # without phylogeny
summary(twobpls2); plot(twobpls2)
# Perform pPCAs
PCAmorph<-phyl.pca(treepls, Morphol, method="lambda", opt="REML")
## [1] 0.4035904 200.6427104
## [1] 0.653023 193.233839
## [1] 0.2494326 203.6558834
## [1] 0.1541578 205.0529098
## [1] 0.09527478 205.74615150
## [1] 0.05888305 206.10832793
## [1] 0.03639173 206.30603125
## [1] 0.02249132 206.41795146
## [1] 0.0139004 206.4831177
## [1] 8.590921e-03 2.065218e+02
## [1] 5.309481e-03 2.065452e+02
## [1] 0.00328144 206.55936689
## [1] 2.028041e-03 2.065680e+02
## [1] 1.253399e-03 2.065734e+02
## [1] 7.746429e-04 2.065767e+02
## [1] 4.787556e-04 2.065787e+02
## [1] 2.958873e-04 2.065799e+02
## [1] 1.828684e-04 2.065807e+02
## [1] 1.130189e-04 2.065812e+02
## [1] 6.984951e-05 2.065815e+02
## [1] 6.984951e-05 2.065815e+02
PCAsound<-phyl.pca(treepls, Sound, method="lambda", opt="REML")
## [1] 0.4035904 208.8991134
## [1] 0.653023 199.919695
## [1] 0.2494326 213.0155644
## [1] 0.1541578 215.1816773
## [1] 0.09527478 216.39513449
## [1] 0.05888305 217.09954352
## [1] 0.03639173 217.51769946
## [1] 0.02249132 217.76955887
## [1] 0.0139004 217.9226938
## [1] 8.590921e-03 2.180164e+02
## [1] 5.309481e-03 2.180739e+02
## [1] 0.00328144 218.10930124
## [1] 2.028041e-03 2.181311e+02
## [1] 1.253399e-03 2.181446e+02
## [1] 7.746429e-04 2.181529e+02
## [1] 4.787556e-04 2.181581e+02
## [1] 2.958873e-04 2.181612e+02
## [1] 1.828684e-04 2.181632e+02
## [1] 1.130189e-04 2.181644e+02
## [1] 6.984951e-05 2.181652e+02
## [1] 6.984951e-05 2.181652e+02
summary(PCAmorph); summary(PCAsound)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 0.05880885 0.04482231 0.02021450 0.01131808 0.007153355
## Proportion of Variance 0.56794037 0.32991774 0.06710322 0.02103597 0.008403048
## Cumulative Proportion 0.56794037 0.89785811 0.96496133 0.98599730 0.994400350
## PC6
## Standard deviation 0.005839448
## Proportion of Variance 0.005599650
## Cumulative Proportion 1.000000000
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 0.03562838 0.02428517 0.01243139 0.003045657
## Proportion of Variance 0.62695850 0.29129216 0.07632838 0.004581508
## Cumulative Proportion 0.62695850 0.91825066 0.99457904 0.999160548
## PC5
## Standard deviation 0.0013036906
## Proportion of Variance 0.0008394516
## Cumulative Proportion 1.0000000000
PC1morph<-PCAmorph$S[,1]
PC1sound<-PCAsound$S[,1]
# PGLS regression between the two PC1s
dataPC1<-cbind.data.frame(names(PC1morph), PC1morph, PC1sound); colnames(dataPC1)[1]<-"Taxon"
dataPCAreg<-comparative.data(treepls, dataPC1, names.col="Taxon")
PCAreg<-pgls(PC1sound~PC1morph, dataPCAreg, lambda="ML")
summary(PCAreg)
##
## Call:
## pgls(formula = PC1sound ~ PC1morph, data = dataPCAreg, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.046548 -0.018683 0.006783 0.034371 0.050127
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.014165
## 95.0% CI : (NA, 0.754)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.0596e-06 6.8984e-02 -0.0001 0.99994
## PC1morph 2.7985e-01 1.5511e-01 1.8042 0.09634 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03289 on 12 degrees of freedom
## Multiple R-squared: 0.2134, Adjusted R-squared: 0.1478
## F-statistic: 3.255 on 1 and 12 DF, p-value: 0.09634
Results are similar to those obtained with our own dataset: the regressions are not significant, probably due to small sample size.
Since phylogenetic signal (l. 140) was high and significant for all
traits except distance TL-labia and IM CSA, we can optimize five out of
seven morphoanatomical traits using contMap.
We had not tested for tracheal elongation ratio, so we need to perform this first.
# Estimating lambda
Tr<-dataS$Tr; names(Tr)<-rownames(dataS)
Tr<-Tr[!is.na(Tr)]
treeTr<-drop.tip(treesyr, setdiff(treesyr$tip.label, names(Tr)))
phylosig(treeTr, Tr, method="lambda", test=TRUE)
##
## Phylogenetic signal lambda : 0.901319
## logL(lambda) : 18.4337
## LR(lambda=0) : 12.2014
## P-value (based on LR test) : 0.000477539
# Evolutionary model selection
for (m in 1:length(models)) {
fit[[m]]=fitContinuous(treeTr, Tr, model=models[m], ncores=2)
}
## Warning in fitContinuous(treeTr, Tr, model = models[m], ncores = 2):
## Parameter estimates appear at bounds:
## a
modTr<-modSel.geiger(fit[[1]],fit[[2]],fit[[3]],fit[[4]],fit[[5]],fit[[6]])
modTr
Phylogenetic signal is high and significant for tracheal elongation ratio. A Brownian Motion model is the best fit.
dataplot=list(); fit=list(); obj=list()
for (i in c(2,4:7)) {
dataplot[[i]]<-dataS[,i]; names(dataplot[[i]])<-rownames(dataS)
dataplot[[i]]<-na.omit(dataplot[[i]])
treeplot<-drop.tip(treesyr, setdiff(treesyr$tip.label, names(dataplot[[i]])))
fit[[i]]<-fastAnc(treeplot, dataplot[[i]], vars=TRUE, CI=TRUE)
obj[[i]]<-setMap(contMap(treeplot, dataplot[[i]]),
colors=rev(brewer.pal(10,"Spectral")))
}
for (i in c(2,4:7)) {
plot(obj[[i]])
title(paste('Ancestral state reconstruction for', colnames(dataS)[i]))
}
dataplot=list(); fit=list(); obj=list()
for (i in c(4:7)) {
dataplot[[i]]<-dataS[,i]/dataS$BM; names(dataplot[[i]])<-rownames(dataS)
dataplot[[i]]<-na.omit(dataplot[[i]])
treeplot<-drop.tip(treesyr, setdiff(treesyr$tip.label, names(dataplot[[i]])))
fit[[i]]<-fastAnc(treeplot, dataplot[[i]], vars=TRUE, CI=TRUE)
obj[[i]]<-setMap(contMap(treeplot, dataplot[[i]]),
colors=rev(brewer.pal(10,"Spectral")))
}
for (i in c(4:7)) {
plot(obj[[i]])
title(paste('Ancestral state reconstruction for', colnames(dataS)[i], 'corrected for body mass'))
}
# With labia length
ggplot(datasyr, aes(x=Family, y=log(Length_labia), fill=Family)) +
geom_boxplot() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
scale_fill_brewer(palette="Dark2")
# With ratio (labia length/body mass)
ggplot(datasyr, aes(x=Family, y=log(Length_labia/Body_mass), fill=Family)) +
geom_boxplot() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
scale_fill_brewer(palette="Dark2")
dataswifts<-datasyr %>% filter(Family=="Apodidae")
datahummingbirds<-datasyr %>% filter(Family=="Trochilidae")
## With labia length
t.test(log(datahummingbirds$Length_labia), log(dataswifts$Length_labia)) # Marginally significant
##
## Welch Two Sample t-test
##
## data: log(datahummingbirds$Length_labia) and log(dataswifts$Length_labia)
## t = 2.5195, df = 6.0742, p-value = 0.04485
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.01949767 1.21182602
## sample estimates:
## mean of x mean of y
## 5.761000 5.145339
## With ratio (labia length/body mass)
t.test(log(datahummingbirds$Length_labia/datahummingbirds$Body_mass),
log(dataswifts$Length_labia/dataswifts$Body_mass)) # Significant
##
## Welch Two Sample t-test
##
## data: log(datahummingbirds$Length_labia/datahummingbirds$Body_mass) and log(dataswifts$Length_labia/dataswifts$Body_mass)
## t = 6.3714, df = 6.0045, p-value = 0.0007001
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.062936 2.388021
## sample estimates:
## mean of x mean of y
## 4.020682 2.295204
familyTER<-c("Caprimulgidae", rep("Apodidae", 3), rep("Trochilidae", 13))
TracheaRatioclean<-na.omit(TracheaRatio)
BMTER<-data$Body_mass[!is.na(data$Tracheal_length)]
TER<-as.data.frame(cbind(familyTER, TracheaRatioclean, BMTER))
colnames(TER)<-c("Family","TER", "Body_mass")
TER[,2]<-as.numeric(TER[,2]); TER[,3]<-as.numeric(TER[,3])
# With TER
ggplot(TER, aes(x=Family, y=log(TER), fill=Family)) +
geom_boxplot() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
scale_fill_brewer(palette="Dark2")
# With ratio (TER/body mass)
ggplot(TER, aes(x=Family, y=log(TER/Body_mass), fill=Family)) +
geom_boxplot() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
scale_fill_brewer(palette="Dark2")
TERswifts<-TER %>% filter(Family=="Apodidae")
TERhummingbirds<-TER %>% filter(Family=="Trochilidae")
## With labia length
t.test(log(TERswifts$TER), log(TERhummingbirds$TER)) # Significant
##
## Welch Two Sample t-test
##
## data: log(TERswifts$TER) and log(TERhummingbirds$TER)
## t = 11.473, df = 13.393, p-value = 2.634e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.4528126 0.6621264
## sample estimates:
## mean of x mean of y
## 4.452572 3.895102
## With ratio (labia length/body mass)
t.test(log(TERhummingbirds$TER/TERhummingbirds$Body_mass), log(TERswifts$TER/TERswifts$Body_mass)) # Not significant
##
## Welch Two Sample t-test
##
## data: log(TERhummingbirds$TER/TERhummingbirds$Body_mass) and log(TERswifts$TER/TERswifts$Body_mass)
## t = 0.93656, df = 3.1209, p-value = 0.4156
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9988831 1.8582229
## sample estimates:
## mean of x mean of y
## 2.135811 1.706141
Gonzalez-Voyer, A. Hardenberg, A. von. 2014. An introduction to
phylogenetic path analysis. In: Garamszegi, L. Z. (Ed.). Modern
phylogenetic comparative methods and their application in evolutionary
biology: concepts and practices. Berlin: Springer, 201–229. https://doi.org/10.1007/978-3-662-43550-2_8
Hardenberg, A. von, Gonzalez‐Voyer, A. 2013. Disentangling
evolutionary cause-effect relationships with phylogenetic confirmatory
path analysis. Evolution 67, 378–387. https://doi.org/10.1111/j.1558-5646.2012.01790.x
Revell, L.J. 2009. Size-correction and principal components for
interspecific comparative studies. Evolution 63, 3258–3268. https://doi.org/10.1111/j.1558-5646.2009.00804.x
Uyeda, J.C., Caetano, D.S., Pennell, M.W. 2015. Comparative analysis
of principal components can be misleading. Systematic Biology 64,
677–689. https://doi.org/10.1093/sysbio/syv019